#include "kernel/chemsys.hpp"
#include "yocto/code/rand.hpp"
#include "yocto/code/ipower.hpp"
#include "yocto/math/kernel/algebra.hpp"

void ChemSys:: randomX() throw()
{
	assert(N>0); assert(M>0);
	assert( X.size() == M );
	for( size_t j=1; j <= M; ++j ) 
		X[j] = 1e-3 + 1.0e-5 * ( 0.5 - alea<Real>() );
}

void ChemSys::compile()
{
	try 
	{
		std::cerr << "---------------- Compiling ChemSys ----------------" << std::endl;

		(size_t &)N = reactions_.size();
		(size_t &)M = library.size();
		
		std::cerr << "-- #Species  = " << M << std::endl;
		std::cerr << "-- #Reactions= " << N << std::endl;
		if( N > M )
			throw exception("*** Too many reactions in ChemSys...");
		if( N <= 0 )
			throw exception("*** No Reactionsin ChemSys!" );
		std::cerr << "-- D.O.F.    = " << M-N << std::endl;
		
		Kproc.free();
		
		
		X.make(M,0);
		Gamma.make(N,0);
		dtGam.make(N,0);
		Nu.make(N,M);
		NuP.make(N,M);
		NuR.make(N,M);
		Phi.make(N,M);
		W.make(N,N);
		iW.make(N,N);
		lss.ensure(N);
		Xi.make(M,N);
		Chi.make(M,M);
		dX.make(M,0);
		
		//----------------------------------------------------------------------
		//-- gather info
		//----------------------------------------------------------------------
		size_t i = 1;
		for( ReactionSet::iterator r = reactions_.begin(); r != reactions_.end(); ++r, ++i )
		{
			Reaction                    &rxn = *r;
			const array<Reaction::Item> &rxn_items  = rxn.items();
			assert( rxn_items.size() > 0 );
			
			//------------------------------------------------------------------
			//-- create function from constant
			//------------------------------------------------------------------
			const Function F( &rxn, & Reaction::getK );
			Kproc.push_back(F);
			
			//------------------------------------------------------------------
			//-- build nu and store integers into NuP and NuR (hack)
			//------------------------------------------------------------------
			for( size_t k=1; k <= rxn_items.size(); ++k )
			{
				const Reaction::Item &it      = rxn_items[k];
				const Species        &species = library[ it.id ];
				const size_t          j       = species.idx;
				assert( it.nu != 0 );
				Nu[i][j] = it.nu;
				if( it.nu > 0 ) 
				{
					*(int *)&NuP[i][j] = it.nu; 
				}
				else
				{
					*(int *)&NuR[i][j] = -it.nu;
				}
			}
		}
		std::cerr << "Nu =" << Nu  << std::endl;
		std::cerr << "---------------- Compiled  ChemSys ----------------" << std::endl << std::endl;

		
	}
	catch(...)
	{
		(size_t &)N = 0;
		(size_t &)M = 0;
		throw;
	}
}

void ChemSys:: computeGammaAndPhi(Real t, bool compute_dtGam)
{
	assert(N>0); assert(M>0);
	
	
	for( size_t i=1; i <= N; ++i)
	{
		//======================================================================
		// outer loop on all reactions
		//======================================================================
		Function   &Kfcn = Kproc[i];
		const Real Ki    = Kfcn(t);
		const Real dotKi = compute_dtGam ? drvs( Kfcn, t, dt ) : 0;
		
		//----------------------------------------------------------------------
		//! compute Gamma and dGamma/dt
		//----------------------------------------------------------------------
		{
			Real       Ri    = 1; //!< initial reactants
			Real       Pi    = 1; //!< initial products
			
			for( size_t j=1; j <= M; ++j )
			{
				const Real   Xj     = X[j];
				const size_t NuR_ij = *(int *)&NuR[i][j];
				const size_t NuP_ij = *(int *)&NuP[i][j];
				
				if( NuR_ij > 0 )
					Ri *= ipower( Xj, NuR_ij);
				
				if( NuP_ij > 0 )
					Pi *= ipower( Xj, NuP_ij);
			}
			Gamma[i] = Ki * Ri - Pi;
			dtGam[i] = dotKi * Ri;
		}
		
		//----------------------------------------------------------------------
		//! compute Phi = dGamma/dX
		//----------------------------------------------------------------------
		{
			for( size_t k=1; k <= M; ++k )
			{
				const Real Xk     = X[k];
				
				Real  Ri = 0;
				Real  Pi = 0;
				const Real NuR_ik = *(int *)&NuR[i][k];
				const Real NuP_ik = *(int *)&NuP[i][k];
				
				if( NuR_ik > 0 )
				{
					Ri = NuR_ik * ipower( Xk, NuR_ik - 1);
					for( size_t j=1; j <= M; ++j )
					{
						if( j != k )
						{
							const Real   Xj     = X[j];
							const size_t NuR_ij = *(int *)&NuR[i][j];							
							
							if( NuR_ij > 0 ) Ri *= ipower( Xj, NuR_ij);
						}
					}
				}
				
				if( NuP_ik > 0 )
				{
					Pi = NuP_ik * ipower( Xk, NuP_ik - 1);
					for( size_t j=1; j <= M; ++j )
					{
						if( j != k )
						{
							const Real   Xj     = X[j];
							const size_t NuP_ij = *(int *)&NuP[i][j];
							
							if( NuP_ij > 0 ) Pi *= ipower( Xj, NuP_ij);
						}
					}
				}
				
				
				Phi[i][k] = Ki * Ri - Pi;
			}
		}
	}
	
	//std::cerr << "X="     << X << std::endl;
	//std::cerr << "Gamma=" << Gamma << std::endl;
	//std::cerr << "dtGam=" << dtGam << std::endl;
	//std::cerr << "Phi  =" << Phi   << std::endl;
	
	
	
}

void ChemSys:: computeXi()
{
	assert(N>0); assert(M>0);
	//----------------------------------------------------------------------
	// Compute W = Phi * Nu'
	//----------------------------------------------------------------------
	algebra<Real>::mul_rtrn(W,Phi,Nu);
	
	
	//std::cerr << "W    =" << W     << std::endl;
	
	//----------------------------------------------------------------------
	// inversion of W by LU decomposition
	//----------------------------------------------------------------------
	if( !lss.LU(W) )
	{
		std::cerr << std::endl;
		std::cerr << "X=" << X << std::endl;
		throw exception("Invalid Chemical Composition!!");
	}
	iW.ld1();
	lss(W,iW);
	//std::cerr << "iW=" << iW << std::endl;
	
	//----------------------------------------------------------------------
	// compute Xi = Nu' * (Phi*Nu') ^{-1}
	//----------------------------------------------------------------------
	algebra<Real>::mul_ltrn(Xi,Nu,iW);
	//std::cerr << "Xi=" << Xi << std::endl;
	
}


void ChemSys:: computeChi()
{
	assert(N>0); assert(M>0);
	//----------------------------------------------------------------------
	// compute Chi = Nu' * (Phi*Nu') ^{-1} * Phi
	//----------------------------------------------------------------------
	algebra<Real>::mul(Chi,Xi,Phi);
	//std::cerr << "Chi=" << Chi << std::endl;
	
}
